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tainties. In addition to the most commonly analyzed data sets for the deep-inelastic scattering of 
charged leptons off nuclei and Drell-Yan di-lepton production, we include also measurements for 
neutrino-nucleus scattering and inclusive pion production in deuteron-gold collisions. The analy- 
sis is performed at next-to-leading order accuracy in perturbative QCD in a general mass variable 
flavor number scheme, adopting a current set of free nucleon parton distribution functions, defined 
accordingly, as reference. The emerging picture is one of consistency, where universal nuclear modi- 
fication factors for each parton flavor reproduce the main features of all data without any significant 
tension among the different sets. We use the Hessian method to estimate the uncertainties of the 
obtained nuclear modification factors and examine critically their range of validity in view of the 
sparse kinematic coverage of the present data. We briefly present several applications of our nuclear 
parton densities in hard nuclear reactions at BNL-RHIC, CERN-LHC, and a future electron-ion 
collider. 
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I. INTRODUCTION AND MOTIVATION 

In spite of the remarkable phenomenological success 
of Quantum Chromodynamics (QCD) as the theory of 
strong interactions, a detailed understanding of the role 
of quark and gluon degrees of freedom in nuclear matter 
is still lacking and subject to ongoing experimental and 
theoretical efforts. In this context, the rather unexpected 
discovery, almost three decades ago, that quarks and glu- 
ons in bound nucleons exhibit non-trivial momentum dis- 
tributions, noticeably different from those measured in 
free or loosely bound nucleons has triggered a long- 
standing quest for more and more precise determinations 
of nuclear parton distribution functions (nPDFs). These 
endeavors have led to increasingly accurate and compre- 
hensive measurements of cross sections involving different 
hard scattering processes and nuclear targets a bet- 
ter theoretical insight into the underlying physics, and a 
more refined framework for analyses of nPDFs (JQ. 

On the one hand, a reliable extraction of nPDFs from 
data is required for a deeper understanding of the mech- 
anisms associated with nuclear binding from a QCD im- 
proved parton model perspective, including a verifica- 
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tion of various proposed nuclear modifications whose phe- 
nomenological details vary from model to model @ , lead- 
ing to a wide spread of expectations. On the other hand, 
nPDFs are a vital input for the theoretical interpretation 
and analyses of a large variety of ongoing and future high 
energy nuclear physics experiments, such as, for instance, 
heavy ion collisions at BNL-RHIC, proton- nucleus col- 
lisions to be performed at CERN-LHC @, Hi , or deep- 
inelastic neutrino-nucleus interactions in long baseline 
neutrino experiments Q . Another important physics ob- 
jective related to nPDFs is to explore and quantify the 
effects of multiple rescatterings and recombinations of 
small momentum fraction gluons, leading to deviations 
fO] from the linear scale evolution usually assumed for 
nPDFs. The transition to the saturation regime is of- 
ten characterized by the saturation scale Q s which de- 
pends on both the relevant momentum fraction x and 
the atomic number A. Important quantitative bench- 
mark tests of saturation phenomena can be performed at 
a future electron-heavy ion collider [ill. [l2j. As a result, 
the kinematic range and the accuracy at which nPDFs 
are known will continue to be a topical issue in many 
areas of high energy nuclear physics. 

In the last few years, significant progress has been 
made in obtaining nPDFs from data. In addition to 
the theoretical improvements nowadays routinely used 
in modern extractions of free proton PDFs, such as the 
consistent implementation of QCD corrections beyond 
the leading order [H and uncertainty estimates [3, Hj], 
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the most recent determinations of nPDFs have also ex- 
tended the types of data sets taken into account, mov- 
in g to wards truly global QCD analyses of nuclear effects 
[5ir fl3l - [l5j . The addition of novel hard probes to the fit 
does not only lead to better constrained sets of nPDFs 
and allows one to study the nuclear modification to the 
different parton species individually, but also tests the 
assumed process independence of nuclear effects. Ver- 
ifying the range of applicability of standard factoriza- 
tion theorems and the universality of nPDFs for hard 
probes in nuclear collisions is of outmost importance as 
formally power suppressed "higher twist" contributions 
can be much enhanced due to the larger density of glu- 
ons in heavy nuclei. 

The deep-inelastic scattering (DIS) of charged leptons 
off nuclear targets not only initiated all studies of nPDFs 
but still provides the best constraints on nuclear modifi- 
cations for quark distributions. Current data comprise a 
wide selection of nuclei from helium to lead, are presented 
as ratios of structure functions for two different nuclei, 
and roughly span the range 0.01<a;<lof momentum 
fractions. Although a separation between quarks and an- 
tiquarks is not possible based on these data alone, DIS at 
medium-to-large x mainly probes valence quarks, while 
data at lower momentum fractions x ~ 0.01 are sensi- 
tive to nuclear modifications of sea quarks. Upon combi- 
nation with available data on Drell-Yan (DY) di-lepton 
production off nuclear targets, a better discrimination 
between valence and sea quarks can be achieved, mainly 
hampered, however, by large experimental uncertainties 
and the limited kinematic coverage. 

DIS and DY data only loosely constrain the nuclear 
modifications to the gluon density because they cover a 
too limited range in the hard energy scale Q such that 
evolution effects are small and at best comparable to the 
current experimental precision. To remedy this situa- 
tion and to further constrain the nuclear gluon density, 
data from BNL-RHIC for inclusive pion production in 
deuteron-gold (dAu) collisions have been included in the 
analysis of nPDFs performed in Ref. [5] . Gluon initiated 
processes are known to be dominant in inclusive hadron 
or jet production at RHIC at not too large transverse mo- 
menta pt [r| , and analogous data for polarized proton- 
proton collisions indeed provide the best constraint on 
the helicity dependent gluon density [17]. Not surpris- 
ingly, the data for dAu collisions at mid rapidity used in 
the fit in Ref. [5J] have a significant impact on their deter- 
mination of the gluon distribution in a gold nucleus. The 
corresponding nuclear modification for gluons at medium 
to large x turned out to be much more pronounced than 
in previous estimates 0, 0] and also much larger than 
those found for all the other partonic species. 

Another promising avenue for significant improve- 
ments in extractions of nPDFs is neutrino induced DIS off 
iron and lead targets available from NuTeV CDHSW 
[l9j , and CHORUS 20] . These data receive their impor- 
tance from their discriminating power between nuclear 
modifications for quarks and antiquarks and have been 



included in a series of analyses in Refs. [13|, |15|. Unex- 
pectedly, the correction factors obtained from neutrino 
scattering data are found to differ significantly both in 
shape and in magnitude from those extracted with the 
more traditional charged lepton probes [l3l [TEj . At vari- 
ance with these results, Ref. [14] confronts the neutrino 
DIS cross sections with nPDFs obtained in Q without 
any refitting and finds no apparent disagreement between 
the nuclear effects obtained from different hard probes. 

The global QCD analysis of nPDFs presented here 
incorporates in a comprehensive way all of the above 
mentioned improvements. The resulting set of nPDFs 
at next-to-leading order (NLO) accuracy supersedes pre- 
vious work presented in [3]. The fitting procedure is 
efficiently performed in Mellin moment space based on 
techniques presented and used in [Tt], l2ll|". W e adopt a 
contemporary set of free nuclcon PDFs 22] defined in 
leading-twist collinear factorization in the MS scheme 
as our reference distribution to quantify modifications 
of PDFs in nuclei. The same general mass variable fla- 
vor number scheme (GM-VFNS) as in [22[ is used in our 
analysis to define charm and bottom quark contributions. 

We utilize the Hessian method [23J to estimate the un- 
certainties of the nuclear modification factors for quarks 
and gluons originating from the experimental errors on 
the fitted data points and examine critically their range 
of validity in view of the sparse kinematic coverage of the 
present data. The resulting eigenvector sets of nPDFs en- 
able one to propagate uncertainties to any desired observ- 
able. We also highlight interesting complications due the 
possibility of having negative parton densities beyond the 
leading order approximation at small values of x in the 
vicinity of typical initial scales of 1 GeV for PDF evolu- 
tion without spoiling the positivity of measured physical 
observables. Such a scenario is realized for our reference 
gluon density at NLO (and beyond) in a free nucleon [22| 
and propagates also to the gluon distribution in nuclei 
obtained in this analysis. 

In addition to the neutral pion production data from 
PHENIX HH used in we include also the charged 
25] and the recently published neutral pion (2(| data 
from the STAR experiment in our fit. Instead of adopt- 
ing only vacuum parton-to-pion fragmentation functions 
(FFs), such as the ones given in Ref. [27j . in our calcu- 
lations, we account for possible medium modifications in 
the formation of the pions by utilizing also a set of nu- 
clear FFs (nFFs) [28| which reproduces the large hadron 
attenuation observed in DIS multiplicities by the HER- 
MES collaboration (29|. Whenever possible, we compare 
measured minimum bias cross sections in dAu collisions 
with our computations at NLO accuracy, rather than uti- 
lizing nuclear modification factors R^Au wnose relation 
to cross sections introduces an additional model depen- 
dence. Regarding the use of neutrino data, we include the 
charged current DIS structure functions F% A and Fg A 
from NuTeV, CDHSW, and CHORUS for iron and lead 
targets [18-20] in our analysis. Mass effects for heavy 
quarks are consistently taken into account using the re- 
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cently obtained expressions of the NLO coefficients [3(| 
in Mellin moment space 31]. They are of particular rel- 
evance for a proper treatment of the strangeness contri- 
bution to charged current DIS which produces a massive 
charm quark in the final state. 

The main features of our new parametrization of 
nPDFs are worth emphasizing already at this point. All 
current data can be described well within conventional 
leading-twist collinear factorization at NLO accuracy by 
a universal set of nPDFs. There are no indications yet 
for the onset of non-linear effects in the scale evolution of 
nPDFs or a breakdown of factorization for hard probes 
involving one heavy nuclei. This is not too surprising 
given the limited kinematic coverage of the data, in par- 
ticular, with respect to the momentum fraction x. We 
find neither the unusually large nuclear modifications of 
the gluon distribution at medium to large x obtained in 
the analysis of [5| nor any tension or discrepancy be- 
tween charged and neutral current DIS results reported 
in [3, EH ■ These differences with previous analyses are 
perhaps a good measure of some, usually disregarded un- 
certainties inherent to global QCD fits such as the ap- 
plied data selection criteria, the different flexibility of 
parameterizations of nuclear modifications, the way of 
propagating experimental uncertainties, or the neglect of 
certain theoretical ambiguities. All these issues need to 
be inspected more closely in the future but likely require 
more precise data and further advances in theory to be 
resolved. 

The remainder of the paper is organized as follows: in 
the next Section we briefly review the general framework 
for a global QCD analysis of nPDFs, establish our con- 
ventions, and describe the strategy of how to parametrize 
nuclear modifications of PDFs in nuclei. In Sec. Ill we 
proceed with a detailed discussion and presentation of the 
results of our analysis. We assess and critically examine 
the uncertainties of nPDFs in Sec. IIIIEI Some expecta- 
tions for future hard probes of nPDFs such as prompt 
photon and DY di-lepton production at RHIC and the 
LHC or hadron yields in DIS are presented in Sec. IIVI 
We summarize the main results of our analysis in Sec. V. 



II. FRAMEWORK 

Throughout this analysis, we make the usual assump- 
tion |3j-|5[ that theoretical expressions for measured cross 
sections da A involving a nucleus A factorize into calcu- 
lable partonic hard scattering cross sections da, identical 
to those used for processes involving free nucleons, and 
appropriate combinations of non-perturbative collinear 
parton densities and fragmentation functions. If applica- 
ble, the latter quantities are subject to nuclear modifica- 
tions and will be denoted as ff and Df' h , respectively. 
Here, i labels the parton flavor, and h represents the 
hadron species produced in the fragmentation process. 
The scale dependence of and Df' h is dictated by the 
proper factorization of collinear mass singularities and, 



hence, will be governed by the same evolution equations 
and kernels as for free nucleons or fragmentation in the 
vacuum. As a result, the entire nuclear modification re- 
sides in the initial conditions for f A and D A ' h at some 
low scale Qo ~ 1 GeV and needs to be parametrized from 
data. 

Applying factorization, the cross sections for DIS, DY, 
and pion production off nuclear beams or targets relevant 
for our global analysis schematically read 

d °ms = J^/i 4 ®^*-** . (!) 

i 

d^y = J^ff^ff^da^ux > (2) 

y 

d^A^x = £ ® // ® da^x ® D A ' n , (3) 

ijk 

respectively, where, for brevity, we have suppressed any 
dependence on kinematic variables, the strong coupling 
a s , renormalization, and factorization scales. We note 
that the DIS cross section da A ls in (JTJ) is usually ex- 
pressed in terms of structure functions L , where I 
denotes either a charged lepton or a neutrino, depend- 
ing on the experimental setup. For the partonic hard 
scattering cross sections da in Eqs. (JJ)-©, the scale evo- 
lution of parton densities and fragmentation functions, 
and the running of a s we consistently use the available 
expressions at NLO accuracy in the MS scheme through- 
out our analysis. We refrain from performing a leading 
order extraction of nPDFs since this leads to a far infe- 
rior description of DY and pion production cross section 
data than in a NLO framework. 

The factorization of all medium related effects into the 
initial conditions for the scale evolution of process inde- 
pendent nPDFs (and, if appropriate, nFFs) is clearly an 
assumption and, apart from its phenomenological suc- 
cess in describing current data, neither proven nor even 
expected to work in general. It provides one, however, 
with a rigorous and testable calculational framework of 
great predictive power, order by order in perturbation 
theory. Global analyses of nPDFs can help to reveal its 
limitations by looking for potential tensions with data. 
Various mechanisms can ultimately lead to a breakdown 
of leading-twist factorization once a regime of dense, sat- 
urated gluons is reached 10]. As a consequence, nuclear 
PDFs are usually applied only to a large class of hard 
probes where a nucleus collides with a lepton, a nucleon, 
or a very light nucleus like the deuteron rather than to 
interactions of two heavy nuclei which create high gluon 
densities. Another reason for neglecting heavy-ion colli- 
sions in nPDF analyses is the complete lack of control of 
the experimentally important impact parameter or cen- 
trality dependence of the probes in a nPDF based frame- 
work. 

The symbol (g> in Eqs. (HJ)-© denotes a convolution 
integral with respect to the relevant momentum fraction. 
To avoid these time consuming integrations in a fit to a 
large body of data, we use the Mellin technique as out- 



4 



lined in Refs. [13, HH which allows one to treat lengthy 
NLO expressions numerically efficient but without resort- 
ing to any approximations. The idea is to represent all 
non perturbative quantities in Eqs. (HJ)- (O) by their rep- 
resentations as Mellin inverses, for instance, 

f t A (x) = ^-.f x~ N tf{N)dN , (4) 

2m J Cm 

where Cat is a suitable contour in the complex N plane 
that has an imaginary part ranging from — oo to +oo and 
that intersects the real axis to the right of the rightmost 
pole of ff(N). Next, after reshuffling integrations in ([!])- 
©, one can compute all quantities, except the desired 
ff~{N) but including the time-consuming da, prior to 
the actual fit and store them in multi-dimensional look- 
up tables in Mellin space [l?], [H|- This technique has 
been successfully exploited in various other global fits of 
parton densities and fragmentation functions Qjl H3, [H, 

At variance with our previous analysis in Ref. ||, 
where the initial nPDFs at scale Qo were related to some 
set of proton distributions through a convolution 

f?(x, Qo) = f A ^-Wf(y, Q )ff (—,Qo) (5) 
Jx N V \ V J 

with appropriately fitted weights W/ 4 , we will work this 
time within a more conventional approach 0, [B[ that de- 
fines the nPDFs for a bound proton in a nucleus A, ff 1 , 
with respect those for a free proton, ff , through a mul- 
tiplicative nuclear modification factor Rf(xN,Qo) as 

ff{x N , Qo) = R?(x N ,Q ) ff(x N ,Q ) . (6) 

Xn resembles the usual DIS scaling variable for free nu- 
cleons, assuming that the momentum of the nucleus pa is 
distributed evenly among its nucleons, i.e., pn — pa/A, 
and has, in principle, support in the interval < xn < A. 
This reflects the fact that a parton in a nucleus may carry 
more than the average nucleon momentum pn . Since the 
ff are restricted to the range < x^ < 1, the nPDFs 
defined through Eq. ([B]) are also constrained to xn < 1- 
Apart from being well suited to Mellin moment space, 
the convolution approach in (J5J) has the advantage to al- 
low for defining nPDFs also beyond xn = 1, see ||. To 
facilitate comparisons to other analyses 0, [H and to em- 
phasize that the results of our fit are not a consequence 
of adopting a different approach, we choose, however, 
the ansatz ([5]) to define our input distributions, which 
has the additional advantage of making the effects of nu- 
clear modifications more transparent than in a convolu- 
tion with a weight function Wf'. 

As the reference PDFs for the free proton, ff, we select 
the latest NLO set from the MSTW global QCD analysis 
22] which is defined in a general mass variable flavor 
number scheme to deal with heavy quark mass effects. 
The nPDFs are then obtained by ^ at an initial scale 
of Qo = 1 GeV by determining the nuclear modification 



factors Rf from data. Their evolution to scales Q > Qo 
follows the prescriptions of the GM-VFNS as specified 
in [22[. This includes the same choices for the initial 
value and the running of the strong coupling and the 
masses and thresholds for the heavy charm and bottom 
quarks. Notice that the medium modified heavy quark 
distributions are generated perturbatively from the gluon 
and light quark flavors, so there is no need to introduce 
any additional free parameters for them. 

Our strategy to parametrize the Rf(xN i Qo) in Eq- © 
is as follows: as in previous analyses 04a], we assume 
isospin invariance and neglect any nuclear effects for the 
deuteron. Both valence quark distributions are assigned 
the same nuclear modification factor Rf(x^, Qo), which 
we parametrize as 

Rf{x,Ql) = e lX a "{l-x)^ x 

(l + e 2 (l-xf*)(l + a v (l-xf*),(7) 

and where we have dropped the subscript N in the parton 
momentum fraction variable for simplicity. The flexible 
functional form in can account for the typical x de- 
pendent pattern of nuclear corrections found in ratios of 
DIS structure functions for different nuclei such as shad- 
owing, anti-shadowing, EMC, and Fermi motion effects. 

We also assume that the light sea quarks and anti- 
quarks share the same correction factor Rf(x,Q^). No 
significant improvement in the quality of the fit to data 
is found by relaxing this assumptions and assigning dif- 
ferent correction factors for each quark flavor. This is 
not too surprising given the limited kinematic coverage 
and precision of the data. We choose another factor 
Rf(x, Qo) to parametrize medium effects for gluons. Ac- 
tually, it turns out that all current data are very well re- 
produced by using nuclear modification factors Rf that 
are not completely independent. An excellent descrip- 
tion of the data is achieved by relating both Rf and Rff 
to Rf specified in Eq. (J7J, allowing only for a different 
normalization and modifications in the low-cc behavior. 
Hence we choose, without any loss in the quality of the 
fit, 

Rf{x,Ql) = Rf( x ,Ql) € -l l + a ^ , (8) 

ei a s + 1 

Rf(x,Ql) = Rf{x,Ql) € JL 1 + a f 9 . (9) 

ei a g + 1 

At large x, where sea quark and gluon densities become 
very small compared to the dominant valence distribu- 
tions, the fit cannot determine extra parameters in their 
nuclear modifications individually. To a first approxi- 
mation it appears to be sensible to assume a common 
large x behavior for all nPDFs. It has to be kept in 
mind, however, that both the resulting EMC effect and 
the Fermi motion for sea quarks and gluons at large x 
are not a result of an experimental constraint but mere 
assumptions which have no impact on the quality of the 
fit. For the use of the Mellin technique outlined above, 
it is important to recall that the N moments of the Rf 
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defined in Eqs. ©-© can be taken analytically, leading 
to appropriate combinations of Eulcr Beta functions. 

To determine the number of actual fit parameters we 
note that the coefficients t\ and £2 in Eq. are fixed 
by charge conservation, i.e., 

f dxf A (x,Q 2 ) = 2 and f dx f£ (x, Q 2 ) = 1 . (10) 
Jo Jo 

This leaves one with nine free parameters per nucleus to 
reproduce all the features of the DIS, DY, and dAu data 
included in the fit, if we further constrain e s and e g to be 
equal, which, again, has no impact on the quality of the 
fit, and fix e s by momentum conservation, 

/ dx y>/^(z,Q 2 ) = i • (11) 

Jo 

The A dependence of the remaining free parameters 
£ € {oiv, a si QW Ah P2, fis, a v , a>si o>g} is parametrized in 
the usual way [3[ as 

£ = 7S + . (12) 

The very mild A dependence found for some of the £'s 
allows us to further reduce the number of additional pa- 
rameters in (fT2")l by setting S ag — 5 as and S ag = 5 Qs , 
leaving a total of 25 free fit parameters. 

The optimum values of these parameters are extracted 
from data by performing a minimization of an effective 
X 2 function that quantifies the goodness of the fit to data 
for a given set of parameters. Given the still sizable ex- 
perimental uncertainties of the data sensitive to nPDFs, 
we choose the simplest \ 2 function, 

i * 

where each experimental result da° xp is compared to its 
corresponding theoretical estimate da th , weighted with 
the uncertainties Aj for each data point. The latter are 
simply estimated by adding statistical and systematic er- 
rors in quadrature. The sum in (|13[) runs over all data 
points i included in the fit, and u)i allows one to give arti- 
ficial weights to different data sets. We refrain from using 
this option and set = 1. In addition, there appears to 
be no need in our fit for introducing relative normaliza- 
tion shifts among different sets of data. We postpone a 
detailed discussion on how we estimate uncertainties of 
our nPDFs to Sec. ITOEl 

III. DISCUSSION OF THE RESULTS 

In this Section we discuss in detail the results obtained 
from our NLO global QCD fit to nuclear scattering data. 
We start with presenting the parameters of the best fit 
and the resulting \ 2 values for each set of data. In the 
following Subsections we discuss the individual probes, 



DIS, DY, neutrino DIS, and dAu collisions, included in 
the fit and show comparisons between data and theory. 
We finish by presenting the obtained nuclear modifica- 
tion factors R A and assessing their uncertainties in Sub- 
sec. HiTeI 

A. Determination of the optimum fit 

The data analyzed comprise the classic EMC [Hj], 
NMC [Mii], and SLAC E139 [13 results for ratios of 
the DIS structure function F A (x,Q 2 ) for various heavy 
nuclei to those for deuterium, lithium, or carbon, see 
Tab. HI We impose a cut Q 2 > 1 GeV 2 on the data to 
ensure that perturbative QCD is applicable, and we are 
in the deep-inelastic regime. We also include the DY di- 
lcpton production data taken in proton-nucleus collisions 
from the E772 [H| and E866 [H collaborations, pre- 
sented as ratios of cross sections for various heavy nuclei 
to those for deuterium and beryllium, respectively. Data 
for single inclusive hadron production in deuteron-gold 
collisions from the PHENIX [H and STAR [H, ^ ex- 
periments are are taken into account for pions at mid ra- 
pidity and px > 2 GeV where NLO QCD provides a good 
description of corresponding pp data (27| . Finally, results 
for neutrino DIS off iron and lead nuclei from the NuTeV 
[3, CDHSW [3, and CHORUS [20[ collaboration are 
included, again after imposing a cut Q 2 > lGeV 2 . The 
total number of 1579 data points considered in our analy- 
sis exceeds those included in our previous fit [3j by almost 
a factor of four, which shows how timely this re-analysis 
is. 

In order to obtain DIS structure functions or parton 
densities for deuterium, needed, for instance, to compute 
pion yields in dAu collisions, we neglect any nuclear ef- 
fects, assume isospin symmetry (u p — d n and d p — u n ), 
and use the free proton PDFs of MSTW [22]. Nuclear 
effects in deuterium were studied in [J] by analyzing data 
on F 2 d /Ff [40] and found be small, 0(1 - 2 %), in partic- 
ular, compared to typical uncertainties of nuclear DIS or 
DY data, which justifies our approach of ignoring them. 
Parton densities in nuclei with A > 2 are constructed 
from the proton densities bound in a nucleus A as de- 
fined in Eq. (fB]), assuming that isospin symmetry also 
holds for bound systems. For instance, the u quark den- 
sity in a nucleus A with Z protons and A — Z neutrons 
at scale fi is given by 

u A (x N ,fi) = jf^{x Nlf i) + ^-j^f£(x N ,fi) , (14) 

and similarly for d A , u A , and d A . Since most data are 
given in terms of ratios of structure functions or cross 
sections, uncertainties in the free proton PDFs in ([B]), 
which can be still substantial for less well determined 
quark flavors or the gluon at small x and low scales Q 
[221 ]. are expected to cancel to a large extent. Neutrino 
induced DIS off nuclei [l8l - l20j is a notable exception since 
the results are presented as absolute structure functions 
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TABLE I: Total and individual \ 2 values for the data sets 
included in the fit. 



measurement 


collaboration 


ref. 


# points 


x 2 




NMC 
E139 


m 

f37] 


17 
18 


18.18 
2.71 


F 2 Ll /F 2 D 


NMC 


f34] 


17 


17.35 


F 2 Li /F 2 D Q 2 dep. 


NMC 


[34] 


179 


197.36 


F 2 Be /F 2 D 


E139 


[37] 


17 


44.17 


F?/F 2 D 


NMC 


[34] 


17 


27.85 




E139 


[37] 


7 


9.66 




EMC 


[33] 


9 


6.41 


F 2 / F 2 Q 2 dep. 


NMC 


[34] 


191 


201.63 


F 2 Al /F 2 D 


E139 


[37] 


17 


13.22 


F? a /Ff 


NMC 


[34] 


16 


18.60 




E139 


[37] 


7 


12.13 


F 2 U 1 F 2 


EMC 


[33] 


19 


18.62 


fP/fP 


E139 


[37] 


23 


34.95 


F 2 g 1 F 2 


E139 


[37] 


7 


9.71 


F 2 Sn /F 2 D 


EMC 


[33] 


8 


16.59 


Fi u /F 2 D 


E139 


[37] 


18 


10.46 


F 2 /F 2 i 


NMC 


[34] 


24 


33.17 


F Ca /F Li 


NMC 


[34] 


24 


25.31 


F 2 /F 2 
F£ l IF§ 
F 2 0o /F 2 ° 


NMC 
NMC 
NMC 


[25] 

[35] 
[35] 


15 
15 
15 


11.76 
6.93 
7.71 


F 2 C 7F 2 C 


NMC 


[35] 


24 


26.09 




NMC 


[35] 


15 


10.38 


F 2 IF 2 


NMC 


[35] 


15 


4.69 


F 2 /F£ Q 2 dep. 


NMC 


[36] 


145 


102.31 


fP/f? 


NMC 


[35] 


15 


9.57 


rpis-t e 
r 2 


NuTeV 


[18] 


78 


109.65 


rni/Fe 
^3 


NuTeV 


[18] 


75 


79.78 


r 2 


CDHSW 


[19] 


120 


108.20 


F3 


CDHSW 


[19] 


133 


90.57 


r 2 


CHORUS 


[20] 


63 


20.42 


^3 


CHORUS 


[20] 


63 


79.58 


j 6' / j D 

dor D y/aa D Y 

j Call D 

da DY /da DY 


E772 
E772 


[38] 
[38] 


9 
9 


9.87 
5.38 


du DY /da DY 


E772 


[38] 


9 


9.77 


d&DY 1 'd^DY 


E772 


[38] 


9 


19.29 


d<TDY/d(TDY 


E866 


[39] 


28 


20.34 


da DY 1 'dooY 


E866 


[39] 


28 


26.07 




PHENIX 


[24] 


20 


27.71 


daT/dall 


STAR 


[26] 


11 


3.92 




STAR 


[25] 


30 


36.63 


Total 






1579 


1544.70 



instead of ratios. In order to account for uncer- 
tainties in the free proton PDFs, we utilize the Hessian 
eigenvector PDF sets of Ref. [22( to estimate the expected 
impact of these PDF variations on the results of our fit to 
^2 3- We add these additional theoretical uncertainties 
in quadrature to the statistical and systematic errors for 

rpvA 
r 2,3 ■ 

The total x 2 f° r the optimum fit was found to be 1544.7 
for 1579 data points and 25 free fit parameters describing 
our nPDFs for quarks and gluons, i.e., a \ 2 per degree 



of freedom very close to unity (% 2 /d.o.f. — 0.994). In 
general, all data sets corresponding to different types 
of observables are adequately reproduced, well within 
the nominal statistical range \ 2 = n ± \/2n with n the 
number of data, cf. Tab. Q] More specifically, the partial 
contribution to \ 2 °f ai l the charged lepton DIS data 
amounts to 897.52 units for 894 data points, for neutrino 
DIS we find 488.20 units compared to 532 data points, 
DY observables amount to 90.72 units for 92 points, and 
pion production in dAu collisions adds another 68.26 
units to x 2 f° r 61 data points. We wish to stress again 
that there is no need to increase the weight uji of any 
particular data set in (|1 3|) to reproduce them well in our 
global analysis. 

The result of the fit suggests that the different data 
sets are complementary in determining the nuclear mod- 
ifications of PDFs for bound protons and that the chosen 
parametrization in Eqs. ©-© and (fT2|) is flexible enough 
to accommodate all the features of the data. We do not 
observe any noticeable tension among the different sets of 
data in the fit. The parameters describing our optimum 
set of nPDFs are listed in Tab .Inland their A dependence 
is illustrated in Fig. [1] We note that the small values for 

for some of the parameters in Tab. HI1 become relevant 
for large A and, hence, are not set to zero. In Tab. lIIII wc 
present for completeness the values for e\^,s,g fixed by 
charge and momentum conservation, Eqs. (fTU|) and (jlip . 
and assuming e s = e g . 



TABLE II: Parameters describing our optimum NLO MS 



nPDFs in Eqs. CO-® and 


l|12[) at the input scale Qo 


= 1 GeV. 


parameter 


7 


A 


S 


a v 


-0.256 


0.252 


-0.017 


a s 


0.001 


-6.89 x 10~ 4 


0.286 


a g 


1.994 


-0.401 


0.286 




-5.564 


5.36 


0.0042 




-59.62 


69.01 


0.0407 




2.099 


-1.878 


-0.436 


a v 


-0.622 


1.302 


-0.062 


a s 


-0.980 


2.33 x 10~ 6 


1.505 


a g 


0.0018 


2.35 x 10~ 4 


1.505 



TABLE III: Values for ei,2,s, ff for our optimum fit in Tab. UTI 
for selected nuclei A as obtained from charge and momentum 
conservation, Eqs. (|10[) and |TT}, and assuming e s = e g . 



A ei e2 e s = e g 

4 0.6612 -0.1033 0.6448 

12 0.7149 -0.1851 0.7147 

27 0.7458 -0.2287 0.7655 

40 0.7596 -0.2487 0.7947 

56 0.7714 -0.2668 0.8239 

197 0.8245 -0.3811 0.9020 

208 0.8280 -0.3912 0.8952 



10 io 2 ^ 10 io 2 A 

FIG. 1: A dependence of the fit parameters according to 
Tab. Ull and Eq. (|12p . Note that ei,2,g, s are fixed by the sum 
rules ((10}, and assuming e 9 = e s . 



10 10 x N 1 

FIG. 2: Data for the DIS structure function ratio F^/F° 
from EMC [13] and NMC [33 ] as a function of momentum 
fraction xn compared with the result of our global fit. 



B. Charge lepton DIS and DY data 

We continue the discussion of the results of our fit with 
a detailed comparison with the available charged lepton 
DIS and DY data, which are the core part of all extrac- 
tions of nPDFs 041 H0. 

Figures M and [3] show the ratios Ff/Fg of the DIS 
structure functions for various nuclei A with respect to 
deuterium from EMC and NMC [H H3 and the E-139 
collaboration (37}, respectively. The solid lines corre- 
spond to the result of the fit at the scale Q 2 of each data 
point, where F^ is computed with our NLO set of nPDFs 
and F® is obtained with the free proton PDFs of |22j. 

Similarly, Figs. [4] and [5] show comparisons with struc- 
ture function ratios using carbon and lithium as reference 
pjil [35l | . The data clearly show the well known regions of 
shadowing, anti-shadowing, and large xn EMC effect for 
Xn 0.05, xjv ~ 0.1, and xn > 0.3, respectively, and are 
in general well reproduced by the fit at NLO accuracy, 
cf. Tab. U for individual \ 2 values. The only exception 
is the low xn behavior of F^ n / ' F® in Fig. [2] however, 
the fit reproduces very well both the ratio F^ jF§ in 
Fig. 03 and F§ '/F? in Fig. [2] We notice, that the strong 
rise of the ratios to values larger than unity as x^ — > 1 
due to Fermi motion is not seen in the analyzed DIS 
data but shows up prominently in low Q 2 experiments, 
see, e.g. [ill ], where target mass corrections are relevant 
though. 

For the imposed cut Q 2 > lGeV 2 , the analyzed 
charged lepton DIS data cover the range xn > 0.01 and 
about a decade in Q 2 for any given bin in xn as is il- 
lustrated in Fig. [5] for the ratio F 2 Sn / Fg [36]. The ob- 
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FIG. 3: The same as in Fig. [2] but now for the SLAC E-139 
data [13] • Note that the multiple points for a given xn have 
different Q 2 values in the range 2 — 10 GeV 2 . 



served, rather moderate Q 2 dependence, best visible for 
the smallest x n , provides some constraint on the nuclear 
modifications of the gluon density through DIS scaling vi- 
olations. Our fit compares well with similar data on the 
Q 2 dependence of the ratios for F^/F? and F 2 C /F 2 D 
[34| . see Tab. [T] Not surprisingly, no substantial differ- 
ences to previous NLO analyses of nPDFs, such as [3j or 
Q , are found in the quality of the description of charged 
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FIG. 5: The same as in Fig. H but now for F 2 A /F 2 " [H and 
Fi/Ff HI from NMC. 



lepton nuclear DIS data. 

Although a distinction of nuclear modifications for 
quarks and antiquarks is not possible based on the 
charged lepton DIS data alone, they provide a valuable 
constraint on nPDFs by mainly probing R v at medium- 
to-large xn and R s at the lowest available momentum 
fractions xn — 0.01. Despite being not too accurate, DY 
di-muon production data [HI, |39| obtained in pA colli- 
sions help to disentangle valence and sea quarks further. 
Again, experimental results are presented as ratios of pA 
cross sections, see Eq. @ , for heavier nuclei and either 
deuterium 38] or beryllium [39] targets and are shown 
in Figs. [7] and |H respectively. 

The E772 and E866 DY data, obtained with an 
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FIG. 6: Data on the Q 2 dependence of ratio Fi n /F^ for 
fixed bins in xn from NMC [3g] compared to the result of our 
global fit. 
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FIG. 7: Ratios of nuclear DY di-muon yields with invariant 
mass M > 4 GeV from E772 [3£| as a function of the parton 
momentum fraction X2 of the nucleus. The solid lines are the 
result of our fit. 



800 GeV proton beam incident on a fixed target and for 
invariant masses M > 4 GeV of the di-muon pair, probe 
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W/Be 




CC DIS off a nucfeon A given, to LO accuracy, by 
F£ A (x N ) - x N [u A + c A + d A + s A }(x N ) , 



FIG. 8: Similar as in Fig.EJbut now for the E866 data [3S| in 
various bins of the invariant mass M of the di-muon pair. 



the range 0.01 < x 2 < 0-2 of momentum fractions x 2 in 
the heavy nuclei. Ratios smaller than unity can be taken 
as an indication of shadowing for sea quark densities in 
nuclei, i.e., R s < 1 at small xn- The relatively large 
scale of the data, set by the invariant mass M of the di- 
muon pair, provides some handle on evolution effects in 
the global fit. 



C. Neutrino induced DIS off nuclear targets 



Charged current (CC) neutrino DIS data [18H20J are 
one of the major additions to our previous analysis in Q 
and are subject to an ongoing discussion [lJ-Qjl about 
their compatibility, or lack thereof, with the neutral cur- 
rent (NC) DIS data discussed in Subsec. IIIIBI A good 
understanding of potential issues with neutrino DIS data 
is also relevant for conventional PDF analyses of free 
protons where they provide a vital constraint on the 
strangeness distribution. 

Neglecting complications due to Cabibbo-Kobayashi- 
Maskawa mixing for the sake of argument, CC DIS data 
draw their relevance for global PDF fits from the differ- 
ent combinations of up-type and down-type quark flavors 
they are sensitive to. With neutrino and antineutrino 
beams one can probe four different structure functions in 



F» A (x N ) ~ x N [u A + c A + d J 



(xn) , 



F 3 * A (x N ) 



[-(u A + c A ) +d A + s A ] (x N ) , 

[u A + c A - (d A + s A )} (x N ) , (15) 



where we have suppressed the scale dependence. Assum- 
ing, as usual, that isospin symmetry holds to a good 
approximation for bound protons and neutrons, the u A 
density in a nucleus A in (jT5j) is given by Eq. (ITi|) and 
similarly for d A ,u A , and d . Experiments extract, un- 
der certain assumptions, averaged structure functions 
F 2 ,3 = (F 2 vA + F 2 A )/2 from appropriate linear combi- 
nations of neutrino and antineutrino CC DIS differential 
cross sections [18l420j . corrected for QED radiative cor- 
rections. As can be easily inferred from (|T5|) . F 2 is pro- 
portional to the total singlet combination of quarks and 
antiquarks and hence sensitive to both valence and sea 
quarks depending on the value of xn. Since the kine- 
matic coverage of CC and NC F 2 data overlaps to some 
extent, any significant tension between these two mea- 
surements should show up prominently in a global QCD 
analysis. Any different nature of interactions of photons 
and charged weak bosons with nuclear matter shall re- 
sult in different, non-universal sets of nPDFs. Likewise, 
the averaged CC structure function F3 mainly probes the 
valence combination u A + d A , which is already well con- 
strained at sufficiently large xn by the NC data discussed 
in Subsec. ImBl By combining CC F 3 and NC F 2 data 
one can arrive at a much improved valence and sea quark 
separation in the entire xn region where data overlap. 

As it turns out, the CC data for the averaged struc- 
ture function F 2 are remarkably well reproduced within 
the experimental uncertainties by our fit, both in shape 
and in magnitude as is illustrated by Fig. [5J The only no- 
ticeable exception are the CDHSW data [19| at Q 2 values 
below 10 GeV 2 where they exhibit a rather different slope 
than the other data. In fact, in this Q 2 region it appears 
to be impossible to simultaneously fit all data sets equally 
well, suggesting some systematic discrepancy among the 
different neutrino data which needs to be further inves- 
tigated. Data for the averaged structure function xnF 3 
are also well described by our fit as can be seen in Fig.fTUl 
except perhaps for the lowest Q 2 values in some bins at 
intermediate xn where the slopes do not match. Both 
sets of data, F 2 and xnF 3 , feature the typical pattern of 
scaling violations, i.e., they increase and decrease with 
Q 2 for x 7v < 0.2 and xn > 0.2, respectively. As men- 
tioned above, the error bars in Figs. [9] and [10] comprise 
the statistical and systematic uncertainties of the data 
added in quadrature and estimates of theoretical ambi- 
guities from variations of the PDFs of free nucleons which 
are most relevant at small xn and low Q 2 (22J. The lat- 
ter are included because neutrino induced DIS data are 
presented as absolute cross section or structure function 
measurements rather than ratios in the absence of vp or 
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FIG. 9: The averaged CC DIS structure function Fa as a function of Q 2 in various bins of xn for iron fl8l . [l9| and lead (20I] 
targets compared with the result of our NLO fit shown as solid and dashed lines, respectively. 



vD benchmark results. Extractions of the nuclear modi- 
fication factors Rf in ([6]) from CC data are hence more 
sensitive to the assumptions about the PDFs of free nu- 
cleons and their uncertainties. 

At variance with our results, a significant tension be- 
tween CC and NC current nuclear DIS data was reported 
in Refs. [HHHI based on a fit of the several thousand data 
points on differential vA and DA cross sections rather 
than the averaged CC structure functions F2.3 used in 
our analysis. Their result, if true, casts serious doubt on 
the validity of pQCD factorization for processes involving 
bound nucleons as it suggests a different, non universal 
behavior of nuclear corrections Rf in CC neutrino in- 
duced and NC charged lepton DIS. This is illustrated, 
e.g., in Fig. 1 of [lj| where the authors compute the ra- 
tio F£ Fe /F£ D from the NuTeV data with an iron target 
and an estimate of the hypothetical structure function 
F% D for neutrino DIS off deuterium, which differs sig- 
nificantly from the measured F^ e jF^ obtained from 
charged lepton DIS data. 

However, our global QCD analysis of nPDFs, in partic- 
ular, the results presented in Figs. l9l and ITOl and Tab. U 
do not support such a strong conclusion which would 
have far reaching implications not only for extractions 



of nPDFs but for free proton PDFs as well as vA DIS 
data for isoscalar targets are often used to constrain the 
strangeness and anti-strangeness densities. We notice 
that also in Ref. [l4j no apparent disagreement between 
the nuclear modification factors Rf for CC and NC DIS 
data has been found based on a comparison of an existing 
fit [5|, not including CC data, with the same set of data 
points for vA and DA cross sections used in [l3|, UH . 

To further illustrate the consistent picture of nuclear 
modifications emerging from our fit, we also show esti- 
mates of the ratio F% e I F^ on the left hand side of 
Fig. [TTJ The experimental results are obtained with the 
NuTeV data closest to the Q 2 values selected in the plot, 
rescaled by a NLO calculation of F% D adopting our ref- 
erence set free proton PDFs from MSTW [22], including 
heavy quark mass effects, and assuming that nuclear ef- 
fects are negligible for deuterium. These data-based ra- 
tios exhibit a pattern of nuclear modifications which is 
not exactly the typical one but resembles all of the ex- 
pected features in that it has shadowing, anti-shadowing, 
EMC, and Fermi motion effects of similar magnitude, at 
low, intermediate, and large values of a; at, respectively. 
The solid lines are computed in a similar way but now 
using the result of our fit to the NuTeV data. These fit- 
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FIG. 10: Same as in Fig. [9] but now for the averaged CC structure function Xjv-Fb. 



driven ratios are clearly consistent with the data-based 
F% Fe jF% D within the present, rather large uncertainties, 
except for the largest values of xjv, x N = 0.65 and 0.75, 
where the theoretical estimate of F^ is most likely not 
reliable enough and where target mass corrections might 
become relevant [3]. We believe that the way in which 
i*2 i s estimated is the main difference with respect to 
what is shown in Il5| . On the right hand side of 
Fig. [H] we present for comparison the corresponding re- 
sults for the nuclear modifications i 7 ^ 6 /i 7 ^ 15 obtained 
from DIS of charged leptons off an iron target. With 
the exception of a slightly more pronounced dip from 
the EMC effect at large xn, the ratios for CC and NC 
DIS are very similar and, unlike Fig. 1 in Ref. [l5l ]. no 
significant tension is observed. A moderate difference be- 
tween the two ratios should be actually expected as they 
probe different combinations of quark densities. How- 
ever, a flexible enough paramctrization of nuclear effects 
Rf can accommodate all sets of data equally well. 

We close the discussion on CC DIS by noticing that the 
proper treatment of heavy quark mass effects is an impor- 
tant asset of our global analysis. The mass dependence 
is fully accounted for by using the recently obtained ex- 
pressions of the NLO coefficients [3(| in Mellin moment 
space [3l| . These corrections are known to be of particu- 



lar relevance for the strangeness contribution to CC DIS 
which produces a massive charm quark in the final state, 
and they have a positive impact on the quality of the fit 
in terms of \ 2 - 



D. Pion production in dAu collisions 

Data for single inclusive pion production at mid rapid- 
ity and high transverse momentum pt in dAu collisions 
at RHIC are the other major addition to our previous 
analysis [1] . Figure [12] shows the neutral and charged 
pion minimum bias production cross sections per nucleon 
for dAu collisions measured by PHENIX [24] and STAR 
[HI Hlj], normalized to the corresponding yields in pp. 
The ratios are obtained for pions at mid rapidity and pre- 
sented as a function of their pt, which also sets the hard 
scale for perturbative calculations using Eq. ([3]). The 
various theoretical curves shown in Fig. [T^] are explained 
and discussed below. 

We are limited to using minimum bias data as collinear 
nPDFs do not exhibit any information on the distribution 
of partons in the transverse plane needed for computa- 
tions of the impact parameter or centrality dependence 
of heavy-ion cross sections. Comparing ratios of mea- 
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FIG. 11: left: data and theory driven estimates of the ratio 
F% Fe / F% D , see text for details, right: nuclear modifications 
from charged lepton DIS, J ? 2 eFe /J ? 2 eD , as obtained from the fit 
to SLAC E-139 data. 



sured minimum bias dAu and pp cross sections avoids 
model dependent estimates of the average number of bi- 
nary nucleon-nucleon collisions (N co u) in a given cen- 
trality class, see, e.g., Ref. [26| for experimental details. 
Notice that even in the absence of nuclear effects, these 
ratios are not necessarily expected to be unity as they can 
be affected by isospin effects such as the smaller density 
of u quarks in a nucleus than in a free proton due to the 
dilution from neutrons. However, noticeable numerical 
effects are only expected for electromagnetic probes like 
prompt photons [42[ which couple directly to the electric 
charge of the quarks, see Sec. IIVI 

In general, results from dAu collisions are significantly 
less straightforward to interpret in terms of nuclear mod- 
ification factors Rf than DIS data. Each value of pt 
samples different fractions of the contributing partonic 
hard scattering processes, integrated over a large range 
of momentum fractions xjv> and convoluted with infor- 
mation on the PDFs of the deuterium. Furthermore, 
since px sets the magnitude for the factorization scale in 
Q, the ratios in Figure [T2l not only reflect the amount of 
nuclear modifications but also their energy scale depen- 
dence. The presence of hadronic probes such as pions, 
charmed mesons, or jets in the final-state, all originating 
from strongly interacting quarks or gluons, leads to ad- 
ditional complications. Apart from the nuclear effects on 
parton densities, accounted for by the nPDFs, the cross 
sections are in principle also sensitive to medium induced 
modifications in the hadronization process. In case of in- 
clusive hadron production and assuming factorizability 
of such medium effects, they should be absorbed into ef- 
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FIG. 12: Ratios of the dAu and pp pion production cross sec- 
tions per nucleon at mid rapidity from the PHENIX [2~i| and 
STAR [25],[2(| collaborations compared to the result of our fit 
using modified [2^ (solid lines) or vacuum [27| (dashed lines) 
FFs. Also shown are calculations using the sets of nPDFs 
from 0] and dotted and dot-dashed lines, repectively. 



fective nuclear parton-to-hadron fragmentation functions 
(nFFs), denoted as D£' 77 in Eq. ©. Modificat ions of 
hadron yields have been found to be quite significant, 
with attenuation effects of up to 50% for heavy nuclei, 
for hadron multiplicities in nuclear DIS by the HERMES 
collaboration 1291 . They have been parametrized in terms 
of nFFs in Ref. [28[ in a NLO anal ysis assuming factoriza- 
tion and based on the HERMES (H DIS and the RHIC 
dAu pion production [241426T ] data. 

The solid lines in Fig. [12] represent the result of our 
best fit of nPDFs using the nFFs of Ref. Q in the cal- 
culation of the dAu cross section in ([3]) and the stan- 
dard vacuum FFs from the global analysis of DSS [27j 
to estimate the pp yields. The fit follows well the rise 
and fall of the ratio at small and high pt, respectively, 
but falls somewhat short in reproducing the enhancement 
found at medium px- Owing to the large experimental 
uncertainties, the x 2 f° r this subset of data is neverthe- 
less very good, Xd.Au/ '#data = 1.12 and slightly better 
than the outcome of an otherwise similar fit using vac- 
uum FFs (dashed lines) where XdAw/^data = 1-37 an d 
Xtot — 1560.5. As it turns out, the use of either nFFs 
or vacuum FFs in the fit has some impact on the shape 
of the different estimates for the ratios in Fig. [T2] In 
any case, the obtained nPDFs are contingent upon the 
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accuracy of the set of FFs used in the analysis and the 
validity of factorization for pion production in dAu colli- 
sions, making DY di-leptons or prompt photons the much 
cleaner but experimentally more demanding probe for 
nPDFs. 

Data for neutral pion yields in dAu collisions were 
first incorporated in the global analysis by EPS [j| and 
found to provide a vital new constraint on Rf u - At vari- 
ance with our approach, the authors in Q disregard any 
medium modifications in the hadronization and, most im- 
portantly, put a large weight oJdAu = 20 on this subset 
of data in the minimization of the x 2 function (TTB")) to 
maximize its impact. They achieve a good description 
of the data as can be seen from the dot-dashed lines in 
Fig. [12] It is not too surprising that the dAu data drive 
rather pronounced modifications of the gluon nPDF in 
their fit, see Fig. 3 in [5], not found in our analysis which 
uses uj dAu = 1. 

In order to better understand the correlation between 
the dAu data and a potentially sizable modification of 
the gluon density, it is instructive to estimate first the 
mean momentum fraction (xn) probed at a given pt- 
The standard way of obtaining (xjy) is to evaluate the 
convolutions in ([3]) with an additional factor on xn in 
the integrand and then divide by the cross section itself, 
see, e.g. Eq. (4) in [43j]. Typically, the (xjv) for pion 
production at mid rapidity at RHIC rises from about 0.05 
at pt — 1 GeV to around 0.3 at px — 15 GeV. Therefore, 
for pt > 1 GeV the cross section is mainly sensitive to 
nPDFs in the anti-shadowing region and at larger pt to 
the suppression due to the EMC effect. 

A better description of the data at intermediate pt in 
where we fall somewhat short, is achieved by an ex- 
traordinarily large enhancement (anti-shadowing) of glu- 
ons at x jv ~ 0.1, where quarks are well constrained by 
DIS and DY data. This is induced by the large weight 
WdAu = 20. In our fit, some part of the enhancement 
in this pt region is provided by a slightly larger gluon- 
to-p ion FFs in a nuclear medium, i.e., D^'^/Dg > 1 
[28l | . At larger values of pt , the region of the EMC effect 
comes into play, where again the quarks are already well 
constrained, and the drop of the data is reproduced by 
an R^ u much less the unity at xjq ~ 0.6. The smaller 
quark-to-pion FFs in a medium, D^' 7r '/D* < 1 [28j. in 
accordance with the large hadron attenuation found by 
HERMES [29], contributes to the behavior of the ratio 
found in our fit. In both regions of a; at, the obtained 
modifications of the gluon nPDF in the EPS fit are much 
more pronounced than the corresponding ones for va- 
lence or sea quarks and not supported by our analysis, 
see Fig. [141 below, based on u<iAu = 1- 

Finally, we refrain from using dAu data obtained at 
forward rapidity by BRAHMS [4J] and STAR 45] which 
show a rapidly increasing suppression of the cross section 
ratios for pt < 2 GeV. Although these data might help to 
further constrain the gluon nPDF down to smaller values 
of x 7v , the theoretical scale uncertainties are very large, 
and the application of pQCD is questionable. Also, little 



is known about FFs and a possible medium modification 
in this kinematic region. 



E. nPDFs and their uncertainties 

Estimating the uncertainties of PDFs and FFs ob- 
tained from global x 2 optimizations has become an im- 
portant objective. The most common and practicable 
approach, the "Hessian method" , explores the uncertain- 
ties associated with the fit through a Taylor expansion 
of x 2 ({£}) around the global minimum Xo({£o}), where 
{£} denotes the set of free parameters of the chosen func- 
tional form at the initial scale Qo and {£o} their values 
for the optimum fit. Keeping only the leading quadratic 
terms, the increase A% 2 can be written in terms of the 
Hessian matrix 



H _ i ay 

13 2 dy l dy j 



(16) 



as 



AX 2 = X 2 ({0) ~ Xo(Uo}) = E H ^Vj (17) 



where {y} are the deviations of the parameters {£} from 
their best fit values, and the derivatives in Eq. (I16[) are 
taken at the minimum. An improved iterative algorithm 
has been devised in [23[ to evaluate the derivatives in 
(|16p reliably in case of very disparate uncertainties in 
different directions of the multi-dimensional space with 
-ZV par parameters describing PDFs or FFs in global QCD 
analyses. We adopt this improved Hessian method also 
in our studies. 

It is convenient to express the Hessian ify in terms of 
its iVp ar eigenvectors and replace the displacements {y} in 
Eqs. (TT^l) and ([T7|) by a new set of parameters {z} related 
to the eigenvector directions. If properly scaled by the 
corresponding eigenvalues, surfaces of constant x 2 turn 
into hyper-spheres in {z} space, and the distance from 
the minimum is given by 



An 



A X 2 = 



(18) 



Within this eigenvector representation {z}, one can 
straightforwardly construct 2iV par eigenvector basis sets 
of nPDFs which greatly facilitate the propagation of 
nPDF uncertainties to arbitrary observables O [Ijjj]. 
These basis sets {S^} correspond to positive and neg- 
ative displacements along each of the eigenvector direc- 
tions by the amount T — y Ax 2 still tolerated for an 
acceptable global fit. To estimate the error AO on a 
quantity O away from its best fit estimate O({£o}) ^ i s 
only necessary to evaluate O for each of the 2N paT sets 
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[H, i.e., 



AO = - 



Y,[0(St)-G(Sr)] 



1/2 



(19) 



The simplicity of this procedure is the main advantage of 
the Hessian approach when compared to the more robust, 
but computationally more involved and less user-friendly 
method based on Lagrange multipliers [4q |. One must 
keep in mind though that the propagation of PDF un- 
certainties in the Hessian method has been derived under 
the assumption that a first order, linear approximation 
is adequate. Of course, due to the complicated nature 
of a global fit, deviations, also from the simple quadratic 
behavior in Eq. (|17p . are inevitable, and error estimates 
based on the Hessian method are not necessarily always 
accurate. 

To initiate the discussions of the nuclear modification 
factors Rf and their uncertainties obtained from the 
global fit, we display in Fig.[T3]our results for gold at the 
initial scale of Qo = 1 GeV 2 , where we assume a common 
Rf = Rf = Rf for all sea quark flavors and Rf v = Rf^ ; 
see Eqs. © - ©, dHJ), and Tab. HI Corresponding re- 
sults evolved to a higher scale of Q 2 = 10 GeV 2 are shown 
in Fig. [2J where we also compare to our previous fit Q 
and the recent analysis of EPS @. To illustrate the A 
dependence, nuclear modification factors in Fig. [TJ] are 
given not only for gold (A = 197) but also for beryllium 
(A = 9), iron (A = 56), and lead (A = 208). 
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FIG. 13: The obtained NLO nuclear modification factors 
Rf u (xN), defined in Eqs. (0 - ©, for gold at our initial 
scale of Qo = lGeV 2 . The inner and outer shaded bands cor- 
respond to uncertainty estimates based on (|19[) for A^ 2 = 1 
and 30, respectively. 



Experimental uncertainties are propagated to the ob- 
tained Rf using the Hessian method outlined above. Ex- 
cursions of the individual eigenvector directions resulting 
in a A% 2 of 1 and 30 units are tolerated and shown as the 
inner and outer shaded bands, respectively, in Figs. [13] 
and [14] In general, Ax 2 = 1 will seriously underesti- 
mate the uncertainties of PDFs and FFs, and, hence, 
a much larger Ay 2 is usually tolerated for an accept- 
able fit [13, HI, I27ll47| . Sophisticated dynamical criteria, 
see, e.g., Ref. [22], have been devised within the Hes- 
sian method to estimate a suitable A% 2 corresponding 
to a, say, 68% (one sigma) confidence level for a given 
fit, but details vary. Reasons for deviating from the de- 
fault A\ 2 = 1 are manifold and can be mainly related to 
uncertainties inherent to the theoretical framework used 
to describe the data, which are notoriously difficult to 
quantify. Examples are the choice of the factorization 
scale, the functional form used to parametrize the PDFs, 
or unavoidable approximations curtailing the available 
parameter space. Given the still rather limited amount 
and kinematic coverage of data taken on nuclear targets 
and their relatively large uncertainties compared to data 
constraining free proton PDFs, we take A% 2 = 30, corre- 
sponding to an increase in \ 2 of about 2%, for an estimate 
of nPDF uncertainties. In any case, it should be kept 
in mind that uncertainty bands are only meaningful for 
combinations of nPDFs and in kinematic regions which 
are actually constrained by data. Therefore, for nPDFs, 
uncertainty estimates below x n — 0.01 should be taken 
with a grain of salt and merely reflect extrapolations of 
the chosen functional form. In our fit this is most appar- 
ent for Rf at small xjv, where the bands shown in Fig. [TBI 
suggest rather small uncertainties, but only charge and 
momentum conservation in Eqs. (|10p and (jlll) provide 
some limited guidance on the behavior at small xn- 

It is worth noticing that the shape and magnitude of 
the nuclear modifications Rf for a given flavor at some 
arbitrary scale Q > Qo depends both on the input distri- 
butions shown in Fig. rj3] and the chosen set of reference 
PDFs for free protons. Scale evolution imprints different 
nuclear effects on individual quark flavors even, as in our 



case, one starts with Rf 



Rf and Ri 



at 



the initial scale Qo- These differences can be quite siz- 
able for the strange and non strange sea quarks as can be 
inferred from comparing Rf and Rf in Figs, [lj and [U 
Even for the nPDFs valence distributions, which evolve 
independently of the quark singlet and the gluon, minute 
differences can be noticed at Q > Qo due to the different 
x shapes of f* and j v d , resulting in R£ ^ Rf . 

The Rf v in Figs. [13] and [JJ] exhibit the expected 
textbook-like behavior of shadowing, anti-shadowing, 
EMC effect, and Fermi motion. In the region constrained 
by DIS data, xn 0.01, uncertainties are small, even for 
the conservative tolerance criterion of A% 2 = 30. Our 
results also agree well with previous determinations of 
Rf in 0,1^. Within uncertainties, there is no evidence 
for any significant anti-shadowing at x n ~ 0.1 also for 
sea quarks. Again, agreement with other fits [H is 
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good, in particular, for R%. Deviations found in the nu- 
clear modifications for strange and light sea quarks at 
Q > Qo are strongly influenced by the shapes of the cor- 
responding distributions in the unbound proton adopted 
in each of the fits, which differ significantly, in particu- 
lar, for the least well determined density ff 22l [47|. To 
some extent differences with previous fits [3|, Q can be 
also attributed to the more flexible functional form © 
to accommodate neutrino DIS data, as well as their im- 
pact on the fit. Another factor is the strong correlation 
with Fl£ at intermediate to large xn- As sea quarks are 
mainly constrained by DIS data at the lowest available 
Xn and DY di-lepton production, uncertainties bands are 
smallest for 0.01 < xn < 0.1, and results for xn < 0.01 
are solely extrapolations. 



As already mentioned in Sec. IIIIDI when discussing 
dAu data, our extracted nuclear modifications for 
the gluon density are expected to differ significantly from 
those determined by EPS ||. Indeed, as can be seen in 
Fig. 1141 we find a much less pronounced anti-shadowing 
region around xn — 0.1 and EMC effect at large xn 



than in the EPS analysis, mainly driven by the way 
in which the dAu data are analyzed; see discussions in 
Sec. IIII Dl above. Differences with our previous fit Q 
are small, however, despite not incorporating any dAu 
data and defining R^ through a convolution with free 
proton PDFs, see Eq. |J5]). Compared to EPS, our best 
fit has significantly less shadowing at Q 2 — 10GeV 2 in 
the unmeasured small xn region, but our uncertainty 
band clearly underestimates the true uncertainties in this 
regime and is biased by the chosen functional form which 
is optimized to provide a good description of the data. 

It is important to notice that despite finding a com- 
paratively moderate nuclear modification of the gluon 
density at Q 2 = 10 GeV 2 in Fig. [H] shadowing is much 
more significant for xn < 0.05 at lower scales, see, e.g., 
Fig. [151 Compared to R^ v , Rf, and also R£, exhibit a 
much more rapid scale evolution at low scales as is il- 
lustrated in Fig. [15] The large suppression of hadron 
yields in dAu collisions at forward rapidities found by 
BRAHMS and STAR @, |H| is essentially sensitive to 
R^ at scales around 1 — 2 GeV 2 and xn ~ 0.001, and 
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FIG. 15: Scale dependence of the valence quark (left) and 
gluon (right) nuclear modification factors for a gold nucleus 
as a function of xn- 
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FIG. 16: Typical nuclear modifications obtained for the per- 
turbatively generated heavy quark parton densities compared 
with those of the gluon and sea quarks for two values of Q 2 . 

can be forced to describe the data without spoiling 
the agreement of the fit with any other data; see also the 
discussion in As mentioned above, we refrain from 
doing so as the applicability of pQCD is not guaranteed 
in this kinematic region and the poor knowledge of FFs 
and possible medium modifications are other obstacles. 
As is also apparent from a comparison of Figs.fTSIandrHl 
the uncertainties on at small xn also rapidly shrink 
under Q 2 evolution, presumably due to gluon radiation 
from quarks at large xn, where they are well constrained. 
This was also observed in the EPS analysis, see Fig. 3 in 
Q , which uses a tolerance level of Ax 2 — 50 for 929 data 
points, compared to our A% 2 = 30 for 1579 data points. 

In order to illustrate also the effective nuclear modi- 
fication for heavy quark flavors, we show in Fig. [15] the 
ratios of the the perturbatively generated charm and bot- 
tom nPDFs in a calcium nucleus and their counterparts 
for a free proton from MSTW 22]. Compared to the 
nuclear modifications of the light sea quarks and gluons, 
which are also displayed in Fig. [111 one finds that R^ a 



and follow closely the xjv-shape of R^ a - Such a 
behavior is expected as heavy quarks are generated by 
gluon splitting without requiring any non-perturbative 
input. The resulting R^ a and R^ a are hence fairly close 
to unity. More pronounced modifications Rg of the glu- 
ons as, for instance, obtained in the EPS analysis [H, 
would imprint larger effects also on R^ and R^. There 
is also an interesting hierarchy in the amount of low Xn 
suppression, which is stronger the lighter the quark is. 

Next, we discuss an important peculiarity of our gluon 
nPDF, not encountered or ignored so far in any of the 
previous analyses EH- While LO PDFs can be 

assigned a physical interpretation as probabilities, at 
NLO and beyond, they become scheme-dependent, non- 
physical quantities. In some recent fits of free proton 
PDFs, including our reference set of MSTW [22|, the 
possibility of negative gluons at small momentum frac- 
tions and scales has been entertained and is actually pre- 
ferred in the best fit of MSTW. The MS evolution ker- 
nels exhibit large order-by-order corrections at small x 
[48j , resulting in rather unstable gluon distributions and 
huge corrections; see Fig. 57 and the detailed discussions 
in Sec. 13 of Ref. [22|. Since quarks tend to rise faster 
due to increasing powers of In x in the splitting function 
P qg at higher orders [48|, gluons compensate for that by 
dimishing at small x and Q 2 . 

As a result, the NLO gluon distribution in the MSTW 
analysis becomes valence-like at Q 2 ~ 2 GeV 2 and neg- 
ative at low x for smaller scales as can be inferred from 
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FIG. 17: Gluon distribution in a free (MSTW Q) and 
bound (gold nucleus) proton at low values of Q 2 . The shaded 
band illustrates the 68% confidence level uncertainties at 
Q 2 = 2 GeV 2 as estimated by MSTW. 
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Fig. [TTl At scales Q 2 > 2GeV 2 , the gluon distribution 
starts to exhibit the well-known strong rise at small x. 
Negative gluons as such are not a problem as long as any 
physical cross section stays positive. The DIS structure 
function Fl is presumably the quantity which is most 
sensitive to its gluon contribution. The corresponding 
gluonic coefficient function also receives large order-by- 
order corrections at small x, which counter the decrease 
of the gluons and stabilize Fl beyond the NLO approx- 
imation (49|. Despite having negative gluons at small x 
and low Q 2 , F L is well behaved for Q 2 >2 GeV 2 at NLO 
and down to even lower Q 2 at NNLO; see Fig. 58 in [22j . 

Since our nPDFs are tied to the free proton PDFs 
of MSTW through the ansatz our NLO nuclear 

gluon density inevitably also turns negative at small xjq 
and Q 2 < 2 GeV 2 as can be seen from the lower panel of 
Fig. [TTl Since R A < 1 for small xn at scale Qo, nPDFs 
are slightly less negative, and we have checked that Fl is 
positive for Q 2 > 2 GeV 2 . In any case, it should be kept 
in mind, that the region x^ % 0.01 is not constrained by 
any experimental result for nPDFs. Also, since g p and 
g A evolve differently with scale, the ratio R g is not well 
defined unless both gluon densities turn positive. For 
Q 2 < 2 GeV 2 and small values of a; at, the ratio R g can 
have nodes for our fit. 

In Fig. [T7J we also illustrate the typical uncertainties 
for the gluon PDF in a free proton in the low Q 2 region 
as estimated by MSTW (shaded band). They turn out to 
be very sizable for x values below a few times 10 -4 . Cor- 
responding uncertainties for the gluon distribution in a 
nucleus are even larger due to the extra uncertainties in- 
troduced by Rf. Since we parametrize the nuclear mod- 
ification factors and estimate their uncertainties with re- 
spect to the best fit of MSTW, a self-consistent propaga- 
tion of uncertainties to the nPDFs is not strictly possible; 
see discussions in Sec. 4 of [5(|. To satisfy constraints from 
momentum and charge conservation, a simultaneous fit of 
free and bound proton PDFs would be necessary which 
does not appear to be feasible at this point given the 
limited experimental information with nuclear targets. 

Finally, for completeness, we have a closer look at the 
actual behavior of the x 2 function (|13[) near its minimum. 
As described above, it is advantageous to work with the 
eigenvector directions {z} of the Hessian matrix, where 
surfaces of constant x 2 are turned into hyper- spheres. 
The contours in Fig. [T5] illustrate the overlap of each 
of the original N pm - = 25 fit parameters {£} listed in 
Tab. HIl with the set of eigenvectors {z}. Eigenvector 1 
corresponds to the largest eigenvalue of the Hessian ma- 
trix, i.e., the direction where x 2 changes most rapidly, 
and number 25 relates to the smallest eigenvalue. One 
can see that several eigenvectors have fairly strong cor- 
relations with only one or a small group of fit parame- 
ters, while others, in particular, those corresponding to 
smaller eigenvalues, overlap with more fit parameters. 
Overall, the result is not quite the ideal case, with a one- 
to-one correspondence between {£} and {z}, thus making 
it difficult to draw conclusions. It basically reflects the 




eigenvector direction 



FIG. 18: Correlations between the fit parameters listed in 
Tab. |TI] and the eigenvector directions of the Hessian matrix 
(see text). 



complicated nature of a global analysis with minimiza- 
tions in a highly correlated, multi-dimensional parame- 
ter space and most likely the still insufficient amount of 
experimental information to clearly pin down differences 
between sea and valence quarks on the one hand and sea 
quarks and gluons on the other. The need to parametrize 
not only the x m shape of the nPDFs but also their A de- 
pendence further complicates the task. 

In Fig. [T5]we investigate the x 2 profiles for each eigen- 
vector direction. We vary one of the parameters {z} at a 
time, keeping all other fixed. Of course, since each eigen- 
vector overlaps, in principle, with all fit parameters, as is 
illustrated in Fig. [T51 the latter are all allowed to change 
in this procedure. The variation is done in such a way 
that a given increase A\ 2 — T 2 is produced, and we com- 
pare the actual behavior (solid lines) for each eigenvector 
direction with the parabolic one (dotted lines) assumed 
in the Hessian approach. One can see from Fig. [18] that 
the quadratic approximation works reasonably well for 
most of the eigenvectors, with only few exceptions, most 
noticeable the profile for the direction corresponding to 
the second largest eigenvalue. Here, T > leads to rather 
tiny changes in \ 2 - From Fig. [T8l one can infer that this 
eigenvector is mainly correlated with the parameters a v , 
a v , and a g . The lack of data at small xn might ex- 
plain to some extent the distorted % 2 profile. Overall, 
we conclude from this exercise that our eigenvector sets 
{S^} for Ax 2 = 30 produce reasonable uncertainty esti- 
mates in the kinematic region constrained by data, i.e., 
for xn > 0.01, with some caveats concerning the flavor 
and the quark versus gluon separation of nuclear modifi- 
cations. 
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FIG. 19: Deviations (solid lines) from the expected parabolic 
behavior (dotted lines) Ax 2 = T 2 for each of the 7V par = 25 
eigenvector directions {z} of the Hessian matrix. 



IV. FUTURE PROBES 

As has become clear from the discussions in the previ- 
ous Section, more data are needed to further our knowl- 
edge of nPDFs to a point where one can address questions 
about possible deviations from linear scale evolution or a 
breakdown of factorization. The biggest obstacle for all 
global analyses of nPDFs is the lack of any DIS collider 
data with heavy ion beams. Measurements of the struc- 
ture functions F 2 and, in particular F L (see [13 for a 
recent study), as well as their scaling violations for var- 
ious nuclei A would constrain the initial conditions for 
nPDFs in a vastly extended range of xjv, similar to the 
one where the partonic structure of free protons is probed 
at present. This would decisively determine also the A 
dependence of nPDFs and, most importantly, challenge 
the currently used theoretical framework in a kinematic 
range where large deviations are expected (lo| . Different 
efforts are currently underway towards a realization of an 
electron-ion collider, see Refs. [ll|; for a status of the 
EIC and LHeC projects, but even in the most optimistic 
scenario it will take at least another decade before first 
data will emerge. 

In the meantime, interesting alternative probes are the 
rapidity dependence of inclusive prompt photon and DY 
lepton pair production in dAu and pPb collisions at RHIC 
and the LHC. In particular, yields at forward rapidities, 
where a large x valence quark in the deuteron (proton) in- 
teracts with a wee, small xn parton in the nucleus, may 
reveal novel aspects of nPDFs. Despite having smaller 
cross sections and being experimentally more challeng- 



ing, these electromagnetic probes have the advantage of 
not exhibiting any sensitivity to nuclear modifications in 
the final state. As we have discussed in Sec. MIDI the 
way how the hadronization process and possible medium 
modifications are modeled can have an impact on the 
obtained nuclear gluon distribution. Prompt photon and 
DY di-lepton production will hence shed light on the con- 
sistency of presently determined nuclear effects. 

As was also mentioned above, and will be seen in our 
results below, the only theoretical complication in analyz- 
ing electromagnetic probes is the presence of potentially 
significant isospin effects due to the direct coupling of 
photons to the electric charge of the quarks. This makes 
one sensitive to, for instance, the smaller density of u 
quarks in a nucleus than in a free proton due to the di- 
lution from neutrons, which was discussed in some detail 
in case of prompt photon production in Ref. [\'2 . Such 
effects need to be taken into account when quantifying 
genuine nuclear modifications for bound protons. 

Prompt photon production has been already advocated 
as a probe of the nuclear gluon density at small xjv in 
Refs. [4J, [5l|. Figure 1201 shows expectations for prompt 
photon yields in dAu and pPb collisions, for central and 
forward (77 — 3) photon rapidities, using our set of nPDFs 
(solid lines) and normalized to the corresponding cross 
section in pp collisions. For comparison, the ratios are 
also computed with the nDS [|| and EPS Q sets of 
nPDFs. All calculations are performed at NLO accu- 
racy [52|. To extract the genuine nuclear modifications, 
the computed ratios should be not compared with unity 
but with the dotted lines which indicate the relevance of 
the isospin effect. The latter curves are obtained with 
free proton PDFs throughout, ignoring any nuclear mod- 
ifications for bound protons, and their deviation from 
unity is solely due to the dilution of the u quark den- 
sity in a neutron-rich nucleus where (A — Z) > Z . To 
facilitate the discussions, Fig. [5T] gives estimates of the 
average momentum fractions (x Pi d) and (xp n .Au) probed 
in the proton (deuteron) and the lead (gold) nucleus for 
the results shown in Fig. [20j The results for (x P} d) and 
) are obta ined in the same way as we have out- 
lined in Sec. HITDl 

At RHIC energies, central rapidity rj — 0, and for 
Pt — 10 GeV, one basically scans the gluon nPDF around 
the anti-shadowing region xm ~ 0.1 and up into the EMC 
effect for larger pt- At r\ — 3, one becomes sensitive to 
smaller momentum fractions but not below xn — 0.01 
already covered by other data. Nevertheless, the pro- 
nounced differences in between the EPS and our fit, 
as illustrated in Fig. [141 lead to characteristic differences 
for daJ Au / daj p , and a measurement at RHIC will cer- 
tainly help to further constrain R^. At the LHC, mo- 
mentum fractions xn down to a few times 10 -3 can be 
accessed with prompt photons produced at forward ra- 
pidities. For r\ — and pr > 20 GeV one mainly probes 
the anti-shadowing peak. Again, any differences between 
the results obtained with the EPS and our set of nPDFs 
in Fig.f2D]are readily explained by the corresponding be- 
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FIG. 20: Expectations for prompt photon production in pPb (left) and dAu (right) collisions at the LHC and RHIC, respectively, 
for central (upper part) and forward (lower part) photon rapidities n using our set of nPDFs (solid lines). Also shown are 
the results obtained with the nDS [|| (dot-dashed lines) and EPS 0] (dashed lines) sets of nPDFs. The dotted lines indicate 
the relevance of the isospin effect (see text). In each case the results are normalized to the corresponding cross section in pp 
collisions calculated with the PDFs from MSTW. 



havior of R A shown in Fig. [14] 

As noticed from our analysis of nPDFs, Drell-Yan pro- 
duction provides an unique tool to disentangle the nu- 
clear effects from valence and sea quark densities. At the 
lowest order in perturbation theory, and keeping only 
the leading u and d quark contributions for the sake of 
simplicity, the (nuclear) cross section ([5} is given by the 
following combination 

da p DY oc e^[u(xi)u A (x 2 ) + u(xi)u A (x 2 )] 

+ e 2 d [d{x 1 )d A (x 2 ) + d(x 1 )d A (x 2 )] • (20) 

Parton distributions are probed at values of x\ i2 which 
depend on the invariant mass M and the rapidity y of the 
dilepton pair (or, equivalently, the intermediate gauge 
boson). Again, at the lowest order the momentum frac- 
tions are given by 

^ 2 = -^e^. (21) 

It follows that at large positive y, where x\ ~ 1 and 
x 2 <C 1, the cross section (l20l) will be dominated by the 
valence distribution of the proton probed at x\ and the 
sea quark u A and d A nuclear modified distributions at 
rather low values of x 2 . The measurement of the cross 



section ratio da^X l^ a ^ provides therefore direct ac- 
cess to the nuclear ratios R A ^{x 2 ). On the other hand, 
at large negative rapidities, distributions are probed at 
ii < 1 and x 2 ~ 1, and Eq. (|2"D|) becomes sensitive to 
the nuclear ratios for the valence distributions instead. 

Figure[^2]shows the expectations for di-lepton pair pro- 
duction in pPb and dAu collisions at the LHC and RHIC, 
respectively, for invariant masses M > 4 GeV and nor- 
malized to the corresponding pp cross sections. The cal- 
culations are performed at NLO accuracy using the code 
from Ref. [53[ and setting the renormalization and factor- 
ization scales as [ip = = M . As in the prompt photon 
analysis, we include the prediction for the ratio using free 
proton PDFs. Isospin effects turn out to be rather small 
for negative rapidities but start to compete with the gen- 
uine nuclear modification at y > (2) for RHIC (LHC) 
energies. It is worth noticing that for RHIC kinematics 
it is possible to cover values of x as low as 10~ 3 at large 
forward rapidity y = 3. In the case of the LHC with 
y/s = 4.4 TeV and at the same rapidity for the di-lepton 
pair, one can explore values of x ~ 5 • 10 ~ 5 , where even 
our present knowledge of the free proton distributions is 
incomplete and will be challenged. It is not unexpected 
then, that it is in this unexplored region where we observe 
the largest differences between the predictions obtained 
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FIG. 21: Estimates of the average momentum fractions (x Pl d) 
and (xp n ,Au) probed in the proton (deuteron) and the lead 
(gold) nucleus for the cross sections shown in Fig. 1201 



with our nPDFs and those of Ref. [5(, which otherwise 
agree due to their similar content of nuclear modification 
in the quark sector, as was already observed in Fig. [TU 
Clearly, measurements of DY cross section ratios at for- 
ward rapidities will further our knowledge of nPDFs. 

Finally, we look into the production of hadrons in 
lepton-nucleus DIS which constitutes an excellent bench- 
mark for different aspects of nuclear effects both in the 
initial and in the final-state. The process is crucially sen- 
sitive to the three main ingredients of a pQCD calcula- 
tion: the effective parton content of the nuclei, the mech- 
anism for partons fragmenting into the detected final- 
state hadron in a nuclear medium, and additional parton 
radiation before and after the interaction with the elec- 
tromagnetic probe. The relevant pQCD framework is 
well known up to 0(a 2 ) [H|, [Hj], and the phenomenologi- 
cal consequences of QCD corrections have been studied in 
detail in [56j. By choosing appropriate kinematical cuts, 
one can enhance different partonic subprocesses which, in 
turn, might be affected differently in a nuclear medium. 

In Ref. [13] the HI collaboration presented a measure- 
ment for neutral pion production in e + p collisions at a 
c.m.s. energy of about y/s — 300 GeV. The 7r°'s were re- 
quired to be produced within a small angle W <E [5°, 25°] 
from the proton beam in the laboratory frame, with an 
energy fraction z v — E^/Ep > 0.01, and 2.5 < pt < 



15 GeV. The data confirmed previous measurements 
which suggested that pQCD predictions at LO accuracy 
underestimate the cross section at low xr [|58lj whereas 
expectations based on BFKL dynamics [59( or on the 
parton content of virtual photons [60j yielded a better 
agreement. The disagreement between the HI data and 
LO estimates based on 0(a s ) cross sections convoluted 
with LO PDFs and FFs could be as large as an order 
of magnitude, depending on the particular kinematical 
region. Since this is much larger than typical higher or- 
der corrections ( U K factor") for such kind of process, it 
was suggested that the data indicate the onset of possi- 
ble non-linear effects in the scale evolution. In Ref. [1^] 
it was pointed out, however, that the particular set of 
kinematic cuts implemented in the HI analysis strongly 
suppresses the LO contributions such that most of the ob- 
served cross section is indeed due to the 7* +g — > g + q + q 
channel, which only opens up at 0(a 2 ), thus explaining 
the large K factor. More specifically, gluon initiated pro- 
cesses in which the pion is produced from a fragmenting 
gluon were found to be dominant. 

Performing similar measurements at a future electron- 
nucleon collider [ll|,[l2j would allow one to probe nPDFs, 
possible medium modifications in the FFs, as well as the 
validity of standard linear scale evolution and collinear 
factorization. In Fig. [231 we show expectations for the 
neutral pion cross section in different bins of Q 2 , obtained 
with different combinations of (n)PDFs and (n)FFs for 
both EIC (left panel) and LHeC (right panel) kinematics. 
In case of an EIC [l2( , we assume collisions of a 30 Gc V 
electron beam with a 100 GeV (per nucleon) gold nucleus 
and similar cuts as in the HI experiment except that the 
transverse momentum of the pion is now allowed to go 
down to 1 GeV. For the LHeC, our results refer to colli- 
sions of 60 GeV electrons with 2.75 TeV per nucleon lead 
ions. Interestingly, for both experiments the predictions 
based on EPS nPDFs @ and ordinary vacuum FFs from 
DSS [27j show a clear enhancement relative to the ex- 
pectations without any nuclear effects based on MSTW 
PDFs, which are largely indistinguishable from the re- 
sults obtained with our new set nPDFs and DSS FFs. 
Using medium modified FFs [28[ instead leads to a sup- 
pression relative to the MSTW(DSS) results. 



V. SUMMARY AND CONCLUSIONS 

We have performed an up-to-date determination of 
parton densities in nuclei [6l| using an extended set of 
data for different observables involving nuclear targets 
and a modern set of parton distributions for free protons 
as reference. The resulting nPDFs are defined at NLO ac- 
curacy in QCD and in a general mass variable flavor num- 
ber scheme. The determination of nPDFs includes error 
estimates obtained within the improved Hessian method 
for A% 2 = 30 and a collection of alternative, eigenvector 
sets of nPDFs that allow one to propagate these uncer- 
tainties, in principle, to any desired observable depending 
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FIG. 22: Rapidity dependence for DY cross section ratios in pPb (left) and dAu (right) collisions at the LHC and RHIC, 
respectively, corresponding to invariant masses M > 4 GeV. The lines represent the expectations for our set of nPDFs (solid), 
EPS [f| (dashed), and the result without considering nuclear effects (dotted). 




FIG. 23: Neutral pion production in DIS in different bins of Q 2 as a function of Bjorken x for four combinations of (n)PDFs 
and FFs. The solid lines represent the result obtained with our best fit of nPDFs and the medium modified FFs of 28]. left: 
Collisions of 30 GeV electrons and 100 GeV per nucleon gold ions at an EIC [12| . The pr of the pion is restricted to be larger 
than 1 GeV. right: Collisions of 60 GeV electrons and 2.75 TeV per nucleon lead ions at the LHeC [llj and pr > 3.5 GeV. 



on these distributions. 

Our results are fully consistent, within uncertainties, 
with a previous determination of nPDFs in Ref. [3] based 
on a much more limited set of data, except for the nuclear 
modifications of the strange quark distribution, mainly 



due to significant changes in the underlying free proton 
reference density. The nuclear modifications for gluons 
are still found to be rather moderate in the entire range 
of momentum fractions, despite including novel exper- 
imental results from dAu collisions. Noticeable devia- 
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tions to Ref. [3j are only found towards small values of 
x and are due to extrapolations outside the kinematic 
region constrained by data. We have also presented nu- 
clear parton densities for charm and bottom quarks, that 
were ignored in the previous analysis. They are generated 
radiatively, i.e., without any additional free parameters, 
from the gluon and light quark distributions in a general 
mass variable flavor number scheme. 

At variance with Refs. [HI, [l5j], we find no conflicting 
patterns of nuclear modifications for neutral and charged 
current deep-inelastic scattering data. We notice, how- 
ever, that the latter set of data have an important impact 
on the shape of the extracted nuclear sea quark densities 
whose modification factors appear to be significantly dif- 
ferent to those found for the valence quarks. 

Compared to the fit in Ref. [5[, which also includes 
some of the available inclusive hadron production data 
in dAu collisions from RHIC, we account also for pos- 
sible nuclear modifications in the hadronization process, 
which arc known to be sizable in multiplicity ratios in 
S1DIS, and refrain from assigning an inflated weight to dAu 
for this subset of data in the fit. The resulting nuclear 
gluon density in the relevant x region constrained by dAu 
data differs considerably from the one obtained in the 
EPS analysis [5]. Compared to the latter fit, which is 
characterized by an anti-shadowing and EMC effect con- 
siderable larger than those for quarks, our gluon den- 
sity exhibits only moderate nuclear corrections. The use 
of standard vacuum rather than modified fragmentation 
functions in our global analysis leads to a marginally 
poorer quality of the fit well inside the tolerated increase 
in x 2 - 

Uncertainties in the nPDFs extraction are still found 
to be rather large, in spite of the inclusion of additional 
data sets, in particular, when compared to the present 
knowledge of PDFs for free protons. As always, the esti- 
mated uncertainty bands depend on the chosen tolerance 
criterion, for which we use A% 2 = 30, and are only trust- 
worthy in the region of momentum fractions constrained 
by data. Extrapolations below xn ~ 0.01 depend mainly 
on the functional form used in the fit and do not reflect 
the true uncertainties. 



Finally, we have presented expectations based on the 
obtained set of nPDFs for some promising hard probes 
comprising prompt photon and forward DY di-lepton 
production at RHIC and the LHC. These measurements 
are expected to further our knowledge of nuclear modifi- 
cations and test their universality. Compared to hadron 
or jet production at RHIC or the LHC, electromagnetic 
probes have the advantage of being independent of possi- 
ble medium modifications of the final state. The biggest 
obstacle in current determinations of nPDFs is, however, 
the complete lack of collider data for DIS off nuclear tar- 
gets, which would constrain nPDFs down to considerably 
lower values of momentum fractions than present fixed- 
target data. It is in this kinematic regime of high gluon 
density where one primarily expects non-linear effects in 
the scale evolution and a breakdown of standard collinear 
factorization. The science case for a future electron-ion 
collider such as an EIC or the LHeC is currently under 
review but first data are at best expected in a decade 
from now. 
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